Multiscale Linear Solvers for Very Large Systems Derived from PDES

We present a novel linear solver that works well for large systems obtained from discretizing PDEs. It is robust and, for the examples we studied, the computational eeort scales linearly with the number of equations. The algorithm is based on a wavelength decomposition that combines conjugate gradient, multi-scaling and iterative splitting methods into a single approach. On the surface, the algorithm is a simple preconditioned conjugate gradient with all the sophistication of the algorithm in the choice of the preconditioning matrix. The preconditioner is a very good approximate inverse of the linear operator. It is constructed from the inverse of the coarse grained linear operator and from smoothing operators that are based on an operator splitting on the ne grid. The coarse graining captures the long wavelength behavior of the inverse operator while the smoothing operator captures the short wavelength behavior. The conjugate gradient iteration accounts for the coupling between long and short wavelengths. The coarse grained operator corresponds to a lower resolution approximation to the PDEs. While the coarse grained inverse is not known explicitly, the algorithm only requires that the preconditioner can be a applied to a vector. The coarse inverse applied to a vector can be obtained as the solution of another preconditioned conjugate gradient solver that applies the same algorithm to the smaller problem. Thus, the method is naturally recursive. The recursion ends when the matrix is suuciently small for a solution to be obtained eeciently with a standard solver. The local feedback provided by the conjugate gradient step at every level makes the algorithm very robust. In spite of the eeort required for the coarse inverse, the algorithm is eecient because the increased quality of the approximate inverse greatly reduces the number of times the preconditioner needs to be evaluated. A feature of the algorithm is that the transition between coarse grids is determined dynamically by the accuracy requirement of the conjugate gradient solver at each level. Typically, later iterations on the ner scales need fewer iterations on the coarser scales and the computational eeort is proportional to N rather than N log N, where N is the number of equations. We have tested our solver on the porous ow equation. On a workstation we have solved problems on grids ranging in dimension over 3 orders of magnitude, from 10 3 to 10 6 , and found that the linear scaling holds. The algorithm works …